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Abstract.  We  study  the  properties  of  a  novel  discontinuous  Petrov  Galerkin  (DPG) 
method  for  acoustic  wave  propagation.  The  method  yields  Hermitian  positive  definite 
matrices  and  has  good  pre- asymptotic  stability  properties.  Numerically,  we  find  that  the 
method  exhibits  negligible  phase  errors  (otherwise  known  as  pollution  errors)  even  in  the 
lowest  order  case.  Theoretically,  we  are  able  to  prove  error  estimates  that  explicitly  show 
the  dependencies  with  respect  to  the  wavenumber  oj,  the  mesh  size  h ,  and  the  polynomial 
degree  p.  But  the  current  state  of  the  theory  does  not  fully  explain  the  remarkably  good 
numerical  phase  errors.  Theoretically,  comparisons  are  made  with  several  other  recent 
works  that  gave  wave  number  explicit  estimates.  Numerically,  comparisons  are  made  with 
the  standard  finite  element  method  and  its  recent  modification  for  wave  propagation 
with  clever  quadratures.  The  new  DPG  method  is  designed  following  the  previously 
established  principles  of  optimal  test  functions.  In  addition  to  the  nonstandard  test 
functions,  in  this  work,  we  also  use  a  nonstandard  wave  number  dependent  norm  on 
both  the  test  and  trial  space  to  obtain  our  error  estimates. 

Key  words:  time  harmonic  wave  propagation;  robustness;  phase  error;  dispersion;  high 
frequency;  Petrov  Galerkin 


1.  Introduction 

The  purpose  of  this  paper  is  to  introduce  a  Discontinuous  Petrov  Galerkin  (DPG) 
method  for  the  Helmholtz  equation.  We  analyze  the  method  and  give  error  estimates 
with  constants  whose  dependence  on  the  wavenumber  are  explicitly  shown.  We  also  report 
results  from  many  numerical  experiments  and  numerically  compare  the  performance  of  the 
DPG  method  with  other  methods.  Although  our  theory  predicts  the  same  h  convergence 
rates  for  the  DPG  method  and  the  standard  FEM,  the  numerical  performance  of  the  DPG 
method  is  far  superior  for  high  wave  numbers.  A  striking  numerical  observation  for  which 
we  do  not  have  a  theoretical  explanation,  is  that  the  DPG  method  exhibits  negligible 
phase  errors. 

The  purpose  of  performing  a  wavenumber  explicit  analysis  is  to  track  down  pollution 
errors  and  gain  a  better  understanding  of  how  they  originate.  These  errors  are  well 
recognized  as  posing  severe  challenges  in  numerical  simulation  of  wave  propagation  [2], 
For  many  model  problems,  the  pollution  is  manifested  as  phase  errors  which  typically 
accumulate  in  the  direction  of  wave  propagation  over  the  computational  domain.  Thus 
the  concepts  of  pollution  error,  phase  error,  and  discrete  wave  numbers  are  all  closely 
related.  To  explain  the  pollution  error  in  the  context  of  finite  element  methods,  we 
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follow  [18]:  Given  that  the  exact  solution  u  lies  in  a  space  U  normed  by  ||  ■  ||t/,  and  the 
discrete  solution  Uh  lies  in  an  approximation  subspace  Uh  C  U,  one  observes  that 

ll",~“l'llg  <  CM  inf  IhfMr,  (1) 

\\u\\u  wheuh  \\u\\u 

where  u  is  the  wavenumber,  C(u)  —  Cj  +  C^cu^cu/i)7,  and  h  is  the  element  size.  The 
inhmum  on  the  right  measures  the  relative  best  approximation  error,  which  is  typically 
controlled  when  uh  is  small,  i.e.,  when  enough  elements  per  wavelength  are  used.  However, 
the  ^-dependence  in  C(u),  as  measured  by  /3,  is  a  reflection  of  the  pollution  errors.  For 
most  standard  methods  of  fixed  order,  the  exponent  f3  is  observed  to  be  a  positive  constant. 

In  contrast,  we  are  able  to  show  that  for  the  new  DPG  method,  quasi-optimality  esti¬ 
mate  (1)  holds  with  C(u)  independent  of  u.  To  be  fully  accurate,  let  us  add  that  we  will 
prove  so  for  an  “ideal”  DPG  method.  The  “practical”  DPG  method  is  a  simple  modifi¬ 
cation  of  the  ideal  DPG  method  (both  methods  will  be  introduced  in  Section  2).  Note 
that  the  independence  of  C(u)  with  respect  to  u  does  not  imply,  in  theory,  that  the  phase 
errors  vanish.  This  is  because  the  U -norm  may  still  contain  cn-dependent  terms.  Yet,  by 
performing  a  wavenumber  explicit  analysis,  we  take  the  first  step  towards  understanding 
why  we  do  not  see  phase  errors  in  our  numerical  experiments.  We  fully  analyze  the  ideal 
DPG  method  and  provide  error  bounds  explicit  in  the  wavenumber  u,  the  mesh  size  h, 
and  the  finite  element  polynomial  degree  p. 

By  doing  so,  we  are  also  able  to  compare  our  work  with  a  few  recent  works  which 
also  give  similar  wavenumber  explicit  estimates.  (These  comparisons  are  in  §  3.2.3.)  For 
perspective,  let  us  recall  the  famous  negative  result  of  [2]  on  the  inevitability  of  pollution 
errors.  Specifically,  Babuska  and  Sauter  [2]  worked  in  the  context  of  a  standard  method 
using  a  nine-point  stencil  (comparable  to  the  standard  FEM  with  bilinear  elements).  One 
then  wonders  what  happens  on  more  general  meshes  and  methods.  It  has  long  been 
known  that  high  order  finite  elements  (with  obviously  larger  stencils)  do  reduce  pollution 
errors.  However,  a  precise  statement  of  this  fact  was  only  recently  obtained  by  Melenk 
and  Sauter  [22],  We  will  compare  our  result  with  this  work.  There  have  been  important 
recent  developments  in  DG  methods  for  the  Helmholtz  equation  [12,  13,  16,  23].  Ultraweak 
formulations,  ever  since  the  works  of  [4,  17],  have  shown  great  potential  in  numerical 
solution  of  the  Helmholtz  equation,  especially  those  approximations  based  on  plane  waves. 
Most  modern  DG  methods  are  derived  from  ultraweak  variational  formulations.  It  is  not 
surprising  therefore  that  a  fertile  line  of  attack  has  been  the  use  of  plane  waves  within 
DG  trial  spaces.  New  theoretical  tools  improving  the  understanding  of  such  methods  have 
just  been  developed  in  the  Plane  Wave  DG  (PWDG)  method  [16].  We  will  compare  our 
results  to  theirs.  Another  method  we  will  compare  with  is  an  interior  penalty  method 
with  complex  stabilization  developed  in  [12,  13].  They  also  provide  an  analysis  of  error 
terms  characterizing  explicitly  the  dependencies  in  wavenumber. 

That  our  method  is  a  Petrov-Galerkin  method  distinguishes  it  from  all  the  above 
mentioned  works.  The  guiding  principle  in  designing  Petrov-Galerkin  methods  is  that 
one  chooses  trial  spaces  for  good  approximation  properties,  but  one  designs  test  spaces 
to  obtain  good  stability  properties.  We  have  exemplified  this  principle  in  our  earlier 
works  [7,  8,  9,  25].  In  [25],  we  considered  the  one- dimensional  Helmholtz  equation  and 
obtained  (1)  with  C(u)  independent  of  u.  Unfortunately  however,  we  were  not  able  to 
generalize  the  techniques  there  to  the  multi-dimensional  case.  In  this  paper,  we  will  pro¬ 
vide  a  different  way  of  theoretical  analysis  that  proves  (1)  for  higher  dimensions  with 
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C(ui)  independent  of  c o.  The  new  analytical  technique  is  built  on  the  approach  developed 
in  [6]. 

We  should  note  that  there  are  two  problems  that  usually  cripple  standard  numerical 
methods  for  the  Helmholtz  equation:  (i)  the  growth  of  the  pollution  error  with  frequency 
and  (ii)  the  poor  approximation  of  highly  oscillatory  wave  solutions  by  polynomials.  One 
can  argue  that  better  trial  spaces  are  needed  to  overcome  the  latter.  Recent  works  such 
as  [16]  raise  the  hope  of  better  trial  spaces.  However,  this  is  not  the  subject  of  this 
paper.  We  concentrate  on  overcoming  the  first  difficulty  by  designing  better  test  spaces. 
Nonetheless,  it  is  important  to  note  that  any  new  developments  in  approximation  spaces 
can  be  built  into  our  approach  easily.  Indeed,  our  method  computes  test  spaces  that  pair 
with  any  given  trial  approximation  space. 

Finally,  let  us  remark  on  a  practically  attractive  feature  of  our  DPG  method:  It  yields 
Hermitian  and  positive  definite  linear  systems  although  the  original  Helmholtz  operator 
is  indefinite.  This  is  not  a  surprise  if  one  views  the  DPG  method  as  a  least  squares 
method.  Indeed,  the  DPG  method  is  a  least  squares  method  in  a  nonstandard  inner 
product,  as  clarified  in  our  earlier  papers  (see  e.g.,  [8,  eq.  (2.13)]).  There  have  been  other 
well  known  least  squares  methods,  such  as  FOSLS,  for  solving  the  Helmholtz  equation, 
notably  [20].  They  show  a  wavenumber- independent  stability  result  if  one  stays  in  a 
subspace  sufficiently  away  from  resonant  modes  (although  it  is  not  clear  how  one  may 
manage  this  numerically).  They  also  claim  reduced  pollution,  but  the  number  of  mesh 
points  they  use  is  far  higher  than  what  we  use  to  obtain  similar  accuracy.  More  interesting 
are  the  results  they  give  on  multigrid  solvers  for  the  least  square  systems  using  the  ideas 
of  [3].  We  hope  to  borrow  these  ideas  to  design  efficient  solvers  for  the  DPG  method,  an 
issue  for  future  research. 

The  paper  is  organized  as  follows.  We  first  describe  the  DPG  method,  in  the  abstract, 
as  well  as  for  the  Helmholtz  application,  in  Section  2.  In  Section  3,  we  present  our  main 
result,  namely  the  wave  number  independent  error  estimate  of  Theorem  3.1.  A  number  of 
comparisons  with  other  recent  works  are  also  made  in  this  section.  In  Section  4,  we  prove 
Theorem  3.1.  In  Section  5  we  present  the  results  of  numerical  experiments  illustrating 
the  theoretical  results  and  comparing  the  DPG  method  to  other  standard  methods. 

2.  The  ideal  and  the  practical  DPG  methods 

In  this  section  we  present  the  new  DPG  method.  We  begin  by  summarizing  the  DPG 
framework  developed  in  [7,  8,  9,  25]  in  §  2.1.  We  then  apply  it  to  the  Helmholtz  setting. 
The  method  we  are  able  to  fully  analyze  is  the  method  presented  in  §  2.2.  The  practically 
implemented  method  is  described  in  §  2.3. 

2.1.  The  abstract  setting.  Let  U  (the  “trial”  space)  and  V  (the  “test”  space)  be  vector 
spaces  over  the  complex  field  C,  and  let  &(•,  •)  :  U  x  V  H »  C  be  a  continuous  sesquilinear 
form.  We  assume  that  U  is  a  reflexive  Banach  space  under  the  norm  ||  •  ||[/  and  that  V 
is  a  Hilbert  space  under  an  inner  product  (•,  -)v  with  a  corresponding  norm  ||  ■  1 1 y .  The 
following  assumption  requires  special  mention,  and  we  will  verify  it  for  the  Helmholtz 
application  later. 

Assumption  1  (Injectivity).  We  assume  that 

{w  E  U  :  b(w,v)  =0,  Vu  <G  V}  —  {0}. 


(2) 
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Now,  suppose  we  are  given  a  continuous  conjugate  linear  form  l(v)  on  V  (we  use  stan¬ 
dard  terminology  -  see  e.g.,  [24]).  The  variational  problem  we  wish  to  approximate  is: 


Find  u  E  U  such  that 
b(u,v)  =  l(v ),  Vu  E  V. 

To  describe  the  DPG  method  for  this  abstract  variational  problem,  we  define 

\b(w,v)\ 


I  opt,  V  =  SUP 
w£U 


m\u 


(3) 


(4) 


and  place  one  more  assumption. 


Assumption  2  (Norm  equivalence).  There  are  positive  constants  Cj ,  C2  such  that 

CilMlv  <  IK'Hopty  <  C*2 ||u|| y ,  Vu  G  V.  (5) 

Clearly,  this  assumption  implies  that  ||u||opt  v  is  a  norm  on  V.  This  is  called  the  optimal 
test  space  norm,  for  reasons  explained  in  [25].  Note  that  the  norms  ||u||opty  and  ||u||v  are 
not  equal  in  general. 

The  approximate  solution  of  the  DPG  method  lies  in  a  finite  dimensional  trial  subspace 
Uh  C  U .  The  test  space  that  pairs  with  the  trial  space  Uh  is  a  subspace  14  C  V  that  we 
now  define:  First,  define  the  trial-to-test  operator  T  :  U  i-»  V  by 

(Tu,v)v  =  b(u,v),  Vu  G  V.  (6) 

Then  the  discrete  test  space  is  set  by 

Vh  =  T(Uh).  (7) 


2.2.  The  ideal  DPG  method.  The  DPG  approximation  Uh  G  Uh  satisfies 

b(uh,vh)  =  l(vh)  Vvhe  14,  (8) 

where  14  is  as  defined  in  (7).  This  is  a  Petrov- Galer kin  type  formulation  as  Uh  and  14 
are  not  generally  identical.  Next,  we  note  two  basic  properties  of  this  method. 

The  first  property  is  that  the  stiffness  matrix  of  the  method  is  Hermitian  and  positive 
definite.  Indeed,  if  {e*}  is  a  basis  for  [4,  then  setting  tj  =  Te*,  we  find  that  the  (i,  j)-th 
entry  of  the  stiffness  matrix,  namely  is 

Ujj  6(ej,4)  (tijtj'iv  (4  Bi)y  bfejUi)  73 '  ji. 

Thus  the  matrix  is  Hermitian.  To  see  that  it  is  also  positive  definite,  let  c  be  a  complex 
vector  and  w  =  c\e \  +  C2e2  +  ...  be  a  basis  expansion  of  any  w  E  Uh-  Then,  by  (6), 

c*Bc  —  b(w,Tw)  —  \\Tw\\y.  (9) 

Now,  by  Assumption  1,  it  is  obvious  that  T  is  injective.  By  (9),  c*Bc  =  0  if  and  only 
if  Tw  —  0,  which  hold  if  and  only  if  c  is  the  zero  vector.  As  a  consequence,  we  can  use 
powerful  iterative  solvers  for  positive  definite  systems  (like  the  conjugate  gradient  method) 
to  obtain  the  DPG  solution,  even  though  the  original  Helmholtz  operator  is  indefinite. 

The  second  property  is  a  basic  convergence  result  for  the  abstract  method.  It  is  also 
easy  to  prove  (see  [25,  Theorem  2.1]),  but  we  omit  the  proof  and  simply  state  it  here. 
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Theorem  2.1:  Suppose  Assumptions  1  and  2  hold.  Then  problems  (3)  and  (8)  are  well- 
posed,  and  their  respective  solutions  u  and  Uh  satisfy  the  quasi- optimality  estimate: 

C2 

\\u  -  Uh\\u  <  yf  inf  || U-WhWu. 

Oi  wheuh 

2.3.  The  practical  DPG  method.  In  view  of  the  fact  that  the  test  space  14  in  (7)  is 
determined  by  T,  it  is  natural  to  ask  if  this  is  computationally  feasible.  As  we  shall  see, 
the  saving  grace  of  the  DPG  formulation  is  that  T  is  a  local  operator  and  consequently 
inexpensive  to  approximate.  Yet,  despite  its  locality,  one  must  approximate  the  infinite 
dimensional  T  by  a  local  finite  dimensional  operator  T  in  a  computer  implementation. 
I.e.,  in  place  of  T,  we  use  T  :  U  eA  V  defined  by 

(Tu,v)v  —  b(u,v),  VugY,  (10) 

where  V  is  the  finite  dimensional  subspace  of  V.  We  will  detail  our  specific  choice  of  T 
and  V  for  the  Helmholtz  application  in  Section  5.  This  perturbation  of  the  ideal  DPG 
method  can  be  analyzed  using  the  recently  developed  techniques  in  [15]  and  will  not  be 
discussed  in  this  paper.  Nevertheless,  it  is  easy  to  verify  that  the  stiffness  matrix  of  this 
practical  method  is  also  Hermitian  and  positive-scmidehnite. 


2.4.  Application  to  the  Helmholtz  equation.  We  consider  a  bounded  domain  Q  C 
Mn  (n  >  2)  with  Lipschitz  boundary.  Let  /  €  L2(i?)  and  g  e  f7_1/2(<9i?).  We  consider  the 
time-harmonic  wave  propagation  problem  as  a  first  order  system.  A  physically  “right” 
way  to  do  this  is  via  the  physics  of  acoustical  disturbances  [5].  Linearizing  the  isentropic 
Euler  equations  around  a  hydrostatic  solution  and  assuming  harmonic  time  variations,  we 
obtain 


iuju  +  V0  =  0, 

on  1? 

(1  la) 

iujcj)  +  V  •  u  =  /, 

on 

(lib) 

u  ■  n  -  0  =  g, 

on  <9l?, 

(11c) 

where  u  and  <f  are  velocity  and  pressure  variables,  respectively,  associated  to  the  acoustic 
perturbations  from  equilibrium.  Observe  that  taking  the  divergence  of  (11a)  and  substi¬ 
tuting  u ,  we  recover  the  usual  second  order  form  of  the  Helmholtz  equation: 


— A0  —  u2cf)  =  iujf 

<90  ... 

- — b  tojcp  =  —  tug 
on 


on  i? 
on  d Q. 


Let  Qh  be  a  disjoint  partitioning  of  i?  into  open  elements  K  such  that  1?  =  U KenhK. 
We  multiply  the  first  two  equations  (11a)  and  (lib)  by  test  functions  and  integrate  by 
parts  element-wise  to  obtain  an  ultraweak  DG  variational  formulation.  The  details  of 
the  derivation  are  very  similar  to  the  case  of  the  Poisson  equation  [6],  so  we  omit  them 
and  simply  present  the  DPG  weak  formulation,  after  a  foreword  on  notations.  Let  (•>  ')d 
denote  the  (sesquilinear)  L2(D )  inner  product  on  any  domain  D.  The  notation  {-,£)i/2,dK 
denotes  the  action  of  a  linear  functional  t  in  H~l/2(dK).  For  concise  notation  that  reflects 
the  element  by  element  calculations,  we  use 


(r,s)t2h  ■=  (r’s)^’ 

KeQh 


(w,£)dnh  ’■=  ^  (■ w,£)i/2,8k ■ 

Kef 4 
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Note  that  in  the  latter  definition,  complex  conjugations  are  absent,  so  to  match  conjugate 
linearity  of  other  terms,  we  will  often  use  notations  like  (w,I)gnh  and  {w,t/QQh,  whose 
meanings  are  self-explanatory.  With  these  notations,  the  equations  of  the  method  derived 
from  the  integration  by  parts  can  be  stated  as  follows: 

iu(u,v)nh  -  {<j>,  V  •  v)nh  +  (<f>,v  ■  n)dQh  =  0  W  G  if  (div,  Qh),  (12a) 

iu(4>,rj)s 4  -  (u,  Vrj)nh  +  ( V,Un)dah  =  (f,v)n,  Vrl  e  H\nh).  (12b) 

Above  and  throughout,  all  derivatives  are  taken  element  by  element  unless  otherwise 
mentioned,  and  the  ‘broken’  spaces  are  defined  by 

H (div,  14)  =  {r  :  t\k  G  H{ div,  K),  \/K  G  14}, 

H\Qh )  =  {u  :  v\K  e  H\K ),  VAT  G  14}- 

From  (12),  it  is  clear  that  there  are  four  solution  components,  namely  ‘interior’  variables 
G  L2(Q)N  x  L2(i ?),  and  the  numerical  trace  and  flux  (fin,  0)  which  lies  in  an  affine 
space  Qg,  which  we  now  define.  Let 

Rg  '=  {{z,n)eH(&\,Q)xHl(Q)\(z-n-n)\dn  =  g}, 

Qg  =  {(2n,P>)  '■  3(z,  n)  G  Rg  such  that  (zn,  jl)  =  tr  dnh(z,p)}, 
where  ( zn,fi )  =  tiQQh{z,p)  signifies  that  for  every  mesh  element  K  G  14,  we  have 

4i|afr  =  •  n\ dK  and  fi\dK  =  p\ ok- 

In  the  case  when  g  —  0,  we  simply  denote  R  =  Ro  and  Q  =  Qo .  Note  that  the  boundary 
condition  (11c)  becomes  an  essential  boundary  condition  imposed  in  the  numerical  trace 
space  Qg.  Observe  that  functions  in  Qgi  when  restricted  to  the  boundary  of  a  single 
element  dK,  are  in  H~lR{dK)xH1^2{dK).  As  already  mentioned,  the  terms  involving  the 
numerical  trace  <f>  and  the  numerical  flux  un  in  (12),  are  to  be  interpreted  as  if_1/2(<91L)- 
functional  actions. 

Let  ( zg,  fig )  in  Rg  and  let  (zg,n,  p>g)  be  its  corresponding  trace  in  Qg.  We  look  for  the 
solution,  decomposed  into 

{u  ,  (j),  Un ,  4)  (iC,  P i  din,  ip)  +  i^Zg,  fig,  ^g,m  Pg)  ■ 

The  component  (zg,  fig,  zg_n,  jig),  consisting  of  the  data  and  its  extension,  is  known.  Hence 
we  only  need  to  compute  the  unknown  (w,  ip,  wn,  dp).  Note  that  (wn,  dp)  has  homogeneous 
boundary  conditions,  i.e.,  it  is  in  Q.  We  can  compute  an  approximation  to  the  unknown 
(w,  ip,  wn,  dp)  by  following  the  abstract  program  in  §  2. 1  §  2.2,  with  these  choices  of  spaces 
and  forms: 


b(  (w,(p,wn,<p), 

(v,v)  )  :=iw(w,v)nh  ~  (p,V  -v)Qh  +  {dp,v  -n)dnh 
+  iu(<p,  rj)(2h  -  (w,  Vv)nh  +  ( V ,  wn)doh, 

(13a) 

K 

(v,v)  )  '■=(f,v)n  ~  {iuzg  +  S7fig,v)n  -  ( iufig  +  V  ■  zg,rj)n, 

(13b) 

U  :=L2(Q)n  x  L2(Q)  x  Q, 

(13c) 

V  :=H(div,f2h)  xH\nh). 

(13d) 

The  norm  on  V  is  defined  by 

II(W^)U2V  =  llVry  +  iuv\\2Qh  +  \\iur]  +  V  -v\\2Qh  +  \\v\\l  +  \\v\\2n, 
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where  the  derivatives  are  calculated  element  by  element  as  usual,  while  in  contrast,  the 
norm  on  R  is  defined  using  the  global  distributional  derivatives: 

II  >  u)\\r  =  ll^h  +  iuz  ||^  +  \\iup  +  V  •  z\\20  +  \\z\\2a  +  WnWjj. 

This  in  turn  defines  the  norm  on  Q  by  standard  quotient  topology,  namely 

\\{zn,fi)\\Q  =  inf  {||(z,//)||fl  :  \/{z,/T)  G  R  such  that  trdnh(z ,  fi)  —  (zn,£i)}.  (14) 

The  U- norm  is  then  inherited  from  the  product  topology,  i.e., 

\\(w,<p,wn,0)\\l  =  \H\l  +  Ml  +  \\(wn,<p)\\2Q-  (15) 

Functions  in  Q  are  single  valued  on  element  interfaces  by  definition.  They  couple  unknown 
interior  values  across  the  mesh  elements. 

3.  The  main  result  and  discussion 

Our  main  result  is  a  wavenumber  independent  quasi-optimality  estimate.  In  this  section 
we  state  this  result  (Theorem  3.1).  Its  proof  follows  from  two  results  proved  in  the  next 
section.  The  remaining  larger  part  of  this  section  is  devoted  to  a  discussion  on  convergence 
rates  and  how  the  method  compares  with  a  number  of  recent  works  by  other  authors. 

3.1.  Wavenumber  independent  quasi-optimality.  Our  analysis  is  based  on  the  fol¬ 
lowing  assumption. 

Assumption  3.  If  f>  satisfies 

on  1?  (16a) 

on  d Q  (16b) 

for  some  F  G  L2(l?)  and  G  G  then  there  is  a  C  >  0  (depending  only  on  Q)  and 

an  ojq  >  0  such  that  for  all  oj  >  ujq,  we  have 

II  W||2D  +  ^U\\l  <  C  (||F||2„  +  ||G||i„)  .  (17) 

This  assumption  is  known  to  hold  on  bounded  convex  domains  (see  [21,  Proposi¬ 
tion  8.1.4]).  It  may  hold  more  generally.  In  fact,  a  more  general  assumption  of  milder 
polynomial  growth  in  the  solution  norm  bound  is  assumed  in  [22,  Assumption  4.7].  How¬ 
ever,  to  our  knowledge,  while  Assumption  3  has  been  verified  for  several  cases,  it  is  still  a 
subject  of  active  research  to  verify  the  more  general  forms  of  this  assumption  on  specific 
domains. 

Theorem  3.1:  Suppose  Assumption  3  holds.  Let  U  =  (w,(p,wn,(p)  G  U  be  the  solution 
of  the  variational  problem  associated  with  the  spaces  and  forms  defined  in  (13)  and  let 
Uh  denote  its  DPG  approximation.  Then  there  exist  constants  ojq  >  0  and  C  >  0  such 
that  the  associated  DPG  solution  Uh  G  Uh  satisfies  the  quasi- optimality  estimate 

\\U-Uh\\u<C  inf  ||U-  Wh\\u,  Vcn>u;o. 

Wh&Uh 

Here,  the  constant  C  is  independent  of  the  wavenumber  c o,  the  mesh  f2h,  and  the  approx¬ 
imating  subspace  Uh-  The  norm  ||  •  ||{/  is  as  in  (15). 

Proof.  This  follows  from  Theorem  2.1,  whose  assumptions  -namely  Assumptions  1  and  2- 
are  verified  in  the  next  section  (see  Theorem  4.5  and  Lemma  4.1).  □ 


A  0  +  u>2(f)  =  F 
^  ±  iujcj)  =  G 
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Note  that  so  far,  we  have  assumed  nothing  about  the  mesh  Q}%  or  the  subspace  U),  C  U 
in  the  above  theorem.  In  fact,  the  theorem  applies  to  arbitrary  element  shapes  and  any 
approximating  subspace  Uh  built  on  any  given  Q p.  However,  it  will  be  useful  to  consider 
a  specific  example  of  Uh  obtained  using  a  tetrahedral  mesh  to  facilitate  comparison  with 
other  works.  We  do  this  next. 

3.2.  Tetrahedral  convergence  rates.  Let  us  now  consider  how  to  use  Theorem  3.1  to 
obtain  convergence  rates  when  1?/,  is  a  geometrically  conforming,  shape  regular,  tetrahe¬ 
dral  finite  element  mesh.  Let  Pp(D )  denote  the  set  of  functions  that  are  restrictions  of 
(multivariate)  polynomials  of  degree  at  most  p  on  a  domain  D.  Let  p  >  0  and  let 

Sh,p  =  {w  :  w\K  G  Pp(K)N},  WKp  =  {v  :  v\K  G  Pp(K )}, 

Qh.,p  Ph)  •  l^h,p)  G  Pi  (ShtP  X  I Vh:p) 

such  that  tr anh(zh,p,  fJ,htP)  =  (znA,ph)}. 

The  example  we  want  to  consider  is  the  case  when  the  trial  space  is  set  by 

Uh  Uh,P  x  kb/jj,  x  Qh,p+ 1-  (IS) 

Clearly,  this  is  a  subspace  of  the  space  U  defined  in  (13c).  We  want  to  derive  h  and  p 
convergence  rates  from  Theorem  3.1.  As  usual,  h  denotes  the  maximum  of  the  diameters 
of  all  mesh  elements.  We  begin  with  the  simplest  case. 

3.2.1.  The  lowest  order  method.  We  set  p  =  0  in  (18)  to  get  the  lowest  order  method,  i.e., 
the  interior  variables  are  approximated  by  piecewise  constant  approximations,  while  the 
numerical  fluxes  and  traces  are  piecewise  linear.  We  only  need  to  study  the  rate  at  which 
the  best  approximation  term  in  Theorem  3.1  converges. 

To  this  end,  let  //,.  denote  the  nodal  interpolant  of  the  linear  Lagrange  finite  element, 
i.e.,  for  a  smooth  function  ip,  the  interpolant  R'lp  on  any  K  G  is  the  linear  function 
whose  values  at  the  vertices  of  K  equal  the  values  of  ip  there.  By  an  abuse  of  notation,  we 
use  the  same  notation  for  vector  functions,  i.e.,  the  vector  function  obtained  by  applying 
Ih  to  each  component  of  v  is  denoted  by  Rv  . 

Let  (u,  <p)  solve  the  Helmholtz  system  (11).  We  tacitly  assume  that  this  pair  is  regular 
enough  to  apply  the  interpolant  Ih-  Since  (u ,  0)  satisfies  the  boundary  condition  of  R ,  the 
approximation  pair  (R u,  Z/i0),  by  construction,  is  also  in  R  (as  can  be  seen  by  comparing 
the  values  of  R u  ■  n  and  <p  at  the  Lagrange  nodes  on  each  face  of  K).  Furthermore, 
t*dnh (IhU,h<f>)  is  in  Qh,P+  i,  so 

inf  \\(un,4>)  ~  (vn,h,'4’h)\\Q  <  \\(un,4>)  ~  tTd{2h(Ihu ,  Ih<t>)\\Q- 
Since  tr  dnh(u,<P)  =  (hn,  0),  by  the  definition  of  the  Q-norrn  in  (14),  we  have 

inf  \\(un,4>)  ~  (vn,h,*Ph)\\Q  <  ||(«,0)  ~  (4«,40)IU-  (19) 

This  is  how  we  bound  the  best  approximation  terms  for  the  numerical  fluxes.  The  other 
terms  forming  the  total  best  approximation  error  in  Theorem  3.1  are  easier.  Combining 
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them  with  (19), 

|| (u,0,0,un)  -  {uh,(f)h,<f>h,unA)\\u  <  C'^||(m»  -  (Ihu,Ih(f))\\2R 

+  Nn|  \\u-wh\\2a+  inf  U  -  iphWh)  ■ 

WhtSh,  o  Vh£Wh,  o  / 

By  the  definition  of  the  i?-norm,  this  implies 

||  ('ll,  0,  0,  fin)  ||  U 

<C  (  uj2\\u  -Ihu\\2n  +  uj2\\(p-  Ih(p\\o  (20) 

+  ||V-(«  ~hu)\\l  +  || V(0  -  h4>)\\2n 

+  \\u-n%u\\2a  +  ||</>  -n^m  ), 

where  77°  and  77{0  denote  the  72-orthogonal  projections  into  Sh,o  and  H4,o,  resp.  In  the 
above  two  inequalities  and  throughout  the  paper  we  use  C  to  denote  a  generic  constant 
independent  of  u.  Its  value  may  differ  at  different  occurrences. 

Convergence  rates  can  now  be  concluded  from  (20).  Note  that  we  used  flux  and  trace 
spaces  of  one  higher  order  than  the  interior  trial  variables.  This  means  that  the  middle 
two  terms  on  the  right  hand  side  in  (20)  converges  at  the  same  h- rate  as  the  last  two 
terms.  By  using  standard  estimates  for  the  L2-projection  and  the  nodal  interpolant,  (20) 
implies  that 

C\\(u,<f>,$,un)  -  (uh,(f>h,$h,unth)\\u  <  u2h2\u\2Hl{Q)  +  uj2h2\(f)\2Hl{n) 

+  hRQ)  +  h2\^\2HHQ) 

+  h-\u\ 2Hl(ty  +  h2\(j)fHl(ay 

At  this  stage  it  is  convenient  to  introduce  a  standard  ^-dependent  norm  (see,  e.g.  [16]), 

m2mD  =  Y,u2('~i)WUoy  (21) 

7=o 

Note  that  if  0  is  a  plane  wave  elL0Xe ,  then  all  the  terms  defining  the  norm  scale  with  u 

in  the  same  way,  namely  as  co2s.  Hence,  ||  •  ||s,w,d  is  often  considered  a  natural  norm  to 

use  for  wave  propagation  problems.  We  use  this  norm  to  summarize  the  conclusion  of  the 
above  discussion. 

Corollary  3.2:  The  DPG  solutions  in  the  lowest  order  tetrahedral  case  satisfy 

II u  -uh ||n  +  || 0  0/i || i?  <  Ch  ( 1 1 0 1 1 T  H|2)W)fl),  (22a) 

||(«n,0)  -  (un,h,j>h)\\Q  <  Ch  ( ||0|| 2,a,,i?  T  \\u  ||2)W,fl)  .  (22b) 

Remark  3.1.  Note  that  although  the  solution  component  u  is  in  77(div,  17),  we  used  the 
nodal  771(i?)-interpolant  to  approximate  it.  We  did  so  only  because  this  is  an  easy  way 
to  find  an  approximating  pair  (Ihu ,  7/,0)  that  satisfies  the  Robin  boundary  condition 
of  R.  If  one  can  find  a  more  natural  interpolant  in  R ,  one  may  be  able  to  improve 
the  regularity  requirements  of  the  above  estimate  (e.g.,  replace  ||w|| 2^,12  by  the  more 
appropriate  norm  ||V  ■  u||i,o;,d)- 
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Remark  3.2.  A  typical  solution  of  the  Helmholtz  equation  is  the  plane  wave  0  =  e  lujd'x 
and  u  —  (f)d  for  some  unit  vector  d  giving  the  direction  of  propagation.  Then,  (22)  implies 
that 

\\u  -  uh\\n  +  ||0  -  <t>h\\n  <  Chu2.  (23) 

This  estimate  shows  that  even  if  uh  is  held  constant,  the  errors  increase  with  u  for  fixed 
data  norms. 

Once  we  fix  uh  to  be  a  constant,  we  may  expect  the  relative  best  approximation  error 
to  remain  more  or  less  constant.  Yet  the  discretization  error  may  not.  (Indeed,  for 
the  example  in  Remark  3.2,  as  we  increase  u,  even  if  we  adjust  h  so  that  uh  remains 
constant,  the  discretization  errors  may  grow  with  u.)  To  our  knowledge  there  is  no 
finite  element  method  to  date  that  can  provably  avoid  this  problem  in  multiple  space 
dimensions.  A  manifestation  of  this  error  increase  with  u,  for  the  standard  methods, 
is  via  the  accumulation  of  phase  errors  in  the  direction  of  wave  propagation.  This  was 
expounded  in  [2]  in  the  context  of  the  standard  method  for  the  Helmholtz  equation,  but 
it  also  well  recognized  for  Maxwell  (see  e.g.  [14,  Fig.  6])  and  other  wave  phenomena. 

Surprisingly  however,  for  the  DPG  method,  phase  errors  were  observed  to  be  negligible 
in  all  our  numerical  experiments.  (See  Section  5  for  an  extended  discussion.)  We  are 
currently  unable  to  explain  this  superior  performance  theoretically.  It  is  possible  that  the 
error  bounds  of  Corollary  3.2  are  too  pessimistic. 

3.2.2.  Higher  order  convergence.  Next,  consider  the  case  p  >  1.  Both  h  and  p  convergence 
rates  can  be  derived  from  Theorem  3.1  because  the  constant  in  the  theorem  is  independent 
of  p.  We  will  now  need  a  conforming  p-optimal  H 1  ( j?)-interpolant,  e.g.,  the  one  given 
in  [11]  and  [10,  Theorem  8.1],  that  satisfies 

110  -  nhp4>\\n  +  h||V(0  -  nhpif)\\Q  <  C  ln(p)2  hs+1p-s\^\Hs+Hnp 

for  all  0  in  77s+1(i?)  with  3/2  <s  +  l<p  +  l.  Above,  p  =  max(p,  2).  For  vector  func¬ 
tions  v  .  let  nhpv  denote  the  vector  function  obtained  by  applying  77 hp  to  each  component 
of  v .  By  following  the  construction  of  77/ip  (see  [11])  it  is  easy  to  see  that  the  approxi¬ 
mating  pair  (nhpu ,  nhp(f>)  satisfies  the  boundary  condition  of  R  whenever  the  pair  (u ,  0) 
is  in  R.  Now,  we  can  proceed  as  in  the  lowest  order  case  (cf.  (20))  to  obtain 

|| u  -Uh\\n+  110-04/2+  ||  (0n,0)  -  {un,h,4>h)\\Q 

<C  ^  u\n(p)2hsp~s\u \HS(I2)  +  ujln(p)2hsp~s\4>\Hs{n) 

+  \n(p)'2hsp~s\u \Hs+i(n)  +  hi(p)2hsp~s\(f)\Hs+i({2) 

+  hsp~s\u  |h»(/2)  +  hsp~s |0|h*(/2) 

Overestimating,  we  can  immediately  summarize  an  estimate  for  the  interior  variables 
using  the  norm  defined  in  (21). 

Corollary  3.3:  The  DPG  solution  in  the  higher  order  tetrahedral  case  satisfies 

ln(h)2 

|| U  —  Uh\\n  +  ||0  —  0/i||/2  +  ChS— — —  (||0||s+l,u;,/2  +  ||^  ||s+l,w,/2)  ,  (24) 


for  all  s  =  1,  2, . . .  ,p. 
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Again,  we  emphasize  that  C  is  independent  of  u,  h,  and  p.  Obviously  a  similar  estimate 
can  also  be  stated  for  the  numerical  traces  and  fluxes. 

3.2.3.  Comparison  with  other  recent  works.  Next,  we  want  to  compare  the  above  stated 
convergence  rates  with  a  few  other  recent  works  that  state  error  estimates  explicitly 
showing  the  wavenumber  dependence.  Note  that  our  method  simultaneously  gives  error 
estimates  for  both  the  pressure  (j)  and  velocity  u .  Most  other  methods  give  only  an 
approximation  to  the  primal  variable  <fi.  An  approximation  to  the  velocity  variable  u 
must  then  be  derived  by  numerical  differentiation  (but  this  results  in  a  loss  of  convergence 
order  for  those  methods).  Hence,  we  will  only  compare  error  estimates  for  <fi. 

(A)  The  standard  p  FEM:  New  estimates  for  this  old  method  have  been  derived  recently 
in  [22],  They  show  that  for  domains  with  analytic  boundary,  or  on  convex  polygons,  if 

—  is  small  and  p  >  Clogu,  (25) 

p 

and  additionally  if  the  Helmholtz  solution  operator’s  norm  satisfies  a  polynomial  growth 
assumption  (satisfied  if  Assumption  3  holds)  then  the  solution  <ph  of  the  standard  /> finite 
element  method  satisfies 

u\\<!>  -  <t>h\\n  +  (ph)\\ah  <  C  inf  (u\\(j)  -  ifh\\n  +  ||V(0  -  tph)\\nh) 

with  a  C  independent  of  c o.  This  is  perhaps  the  clearest  precise  statement  available  in 
the  literature  demonstrating  that  pollution  effects  are  removed  in  high  order  p  FEM.  This 
estimate  is  better  than  the  estimate  of  our  Theorem  3.1.  Yet,  the  numerical  performance 
of  our  method  (in  the  case  of  low  p,  as  reported  later)  is  better  than  the  standard  FEM. 
The  main  advantage  in  our  theory  is  that  we  have  no  need  for  condition  (25).  The  DPG 
method  has  better  pre-asymptotic  stability  properties  (e.g.,  we  have  no  need  to  assume 
a  sufficiently  small  h)  and  yields  Hermitian  positive  definite  matrices.  The  growth  of 
conditioning  with  h  of  both  methods  are  similar. 

(B)  The  plane  wave  DG  method  (PWDG):  In  two  space  dimensions,  the  recent  pa¬ 
per  [16,  Theorem  3.14]  analyzes  a  Trefftz  DG  method  using  plane  waves  for  trial  sub¬ 
spaces.  If  plane  waves  in  p'  =  2m  +  1  wave  directions  (sufficiently  separated)  are  used 
with  each  mesh  element,  then  for  sufficiently  large  p',  they  prove  that 

<v||0  -  4>h\\n  <  C(uh )  diam(h?)  hs~l  H\\s+i,u,o  (26) 

holds  for  all  s  <  \{m  +  l)/2],  where  C{uh )  is  an  increasing  function  of  uh.  For  the 
sake  of  comparison  when  uh  is  fixed,  we  may  multiply  both  sides  of  (26)  by  h  so  that 
both  (26)  and  our  estimate  (24)  gives  the  same  /?, -convergence  rate.  If  one  agrees  to  view 
our  polynomial  degree  p  to  be  more  or  less  comparable  to  their  parameter  p',  then  our 
estimate  is  comparable  to  theirs  (with  the  difference  that  we  neither  have  u  dependence 
in  C,  nor  do  we  need  to  assume  large  p).  However,  since  we  work  with  polynomial  spaces, 
we  avoid  the  conditioning  problems  they  faced  due  to  the  use  of  plane  waves. 

We  should  note  however  that  the  numerical  results  reported  in  [16]  are  excellent.  It 
begs  the  question  if  we  could  use  their  same  plane  wave  basis  functions  to  form  the  trial 
space  Uh  in  our  DPG  setting.  Indeed,  this  can  be  done.  Note  that  in  Theorem  3.1  we 
placed  no  assumptions  on  Uh ,  so  the  theorem  applies  verbatim  in  this  setting.  The  only 
theoretical  difficulty  is  in  bounding  the  best  approximation  error  estimate.  While  the  best 
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approximation  estimates  for  the  interior  variables  immediately  follow  from  [16],  bounding 
the  flux  best  approximation  terms  seems  to  require  conforming  plane  wave  approximation 
estimates  not  available  in  [16]. 

(C)  Interior  penalty  DG  method  with  complex  stabilization:  This  method  was  recently 
developed  in  [12,  13].  In  the  lowest  order  case  [12,  Theorem  5.5  and  Eq.(6.6)]  they  prove 
that 

||V(0-<^)|k  <  (l  +  toh)(C1toh  +  C2to3h2) 

for  a  specific  choice  of  their  stabilization  parameters.  Here  C\  and  C2  depend  only  on  the 
load  and  are  independent  of  to  and  h.  We  may  compare  this  to  our  estimate  (23).  If  c oh 
is  held  fixed,  then  the  growth  with  respect  to  to  in  both  the  estimates  are  linear.  In  [13], 
the  higher  order  case,  for  a  general  polynomial  degree  p,  is  considered.  The  best  estimate 
they  have,  after  an  iterative  improvement  [13,  Theorem  5.1],  is 

Lmin(p+l,s) 

\\<j>  —  <j>h\\f2  <  C - - - ||0|| H°(Q),  (27) 


provided 


V 

The  estimate  (27)  is  comparable  to  (24).  A  notable  difference  is  the  absence  of  factors  of 
ln(p)  in  their  estimate:  These  factors  arose  due  to  our  need  to  use  conforming  /^-optimal 
projectors,  a  need  absent  in  traditional  DG  analyses.  Also  note  that  we  have  no  need  for 
assumption  (28). 


4.  Analysis 

In  this  section,  we  verify  the  assumptions  of  Theorem  2.1  for  the  DPG  method  applied 
to  the  Helmholtz  equation.  Throughout  this  section,  we  tacitly  assume  that  12  is  such 
that  Assumption  3  holds  for  Helmholtz  solutions.  We  show  that  Assumption  3  implies 
Assumptions  1  and  2  for  the  Helmholtz  application. 

4.1.  Verification  of  injectivity.  To  prove  the  following  lemma  we  only  use  a  weak 
consequence  of  Assumption  3,  namely  that  (17)  implies  uniqueness  of  solutions  for  the 
Helmholtz  problem  (16). 

Lemma  4.1:  Assumption  1  holds  for  the  DPG  sequilinear  form  defined  in  (13). 

Proof.  Consider  any  (u  ,<j>,un,^>)  G  U  —  L2(i 1)N  x  L2(i 7)  x  Q  satisfying 

iuj(u,v)nh  -  (</>, 'V  ■  v)ah  +  -n)eQh  =  0, 

-  _  (29) 

M<i>>v)nh  -  {U,  Vr?)^  +  (77,  un)dQh  =  0, 

for  all  (77, 77)  G  V  =  H( div,  Qh)  x  771(i7/l).  Testing  with  functions  in  the  subspace  of  V 
consisting  of  globally  infinitely  differentiable  functions,  compactly  supported  on  i?,  we 
find  that 

itou  +  V0  =  0  and  iuxf  +  V  •  u  —  0  (30) 

in  the  sense  of  distributions  on  the  open  set  12.  This  implies  that  <f>  G  H1(f2)  and 
u  G  i7(div,  12),  which  now  allows  us  to  integrate  by  parts  in  (29).  Hence,  for  every 
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(v,  rj)  e  V,  we  obtain  the  equations 


((ft  —  (ft,  v  ■  n)  =0  and  (77,  un  —  u  ■  n)dQh  =  0. 

Thus,  (un,<t>)  =  Ggah(u ,  (ft).  Furthermore,  since  ( un,(ft )  £  Q,  we  satisfy  the  boundary 
condition 

u  ■  n  —  (ft  —  0,  over  d Q.  (31) 

Now,  we  test  equation  (29)  with  the  globally  conforming  rj  =  (ft  e  Hl(Q)  to  get 

MWn  ~  (u,  V(ft)n  +  UWln  =  °- 

Since  V(ft  =  —itou  by  (30),  the  real  part  of  this  equation  implies  that  =  0.  ffence, 

by  (31),  (ft  =  u  ■  n  —  0  on  dFl.  Thus,  (ft  satisfies,  in  a  distributional  sense,  the  Helmholtz 
boundary  value  problem  (16)  with  zero  F  and  G.  Therefore,  by  Assumption  3,  (ft  =  0  in 
i2.  Then,  we  obviously  also  have  u  —  0  in  i2,  (/)\qk  =  <ft\dK  =  0  and  un\dK  =  u  ■  n\gK  =  0 
for  all  K  G  Qh.  □ 

4.2.  Optimal  test  norm.  The  optimal  norm  is  easily  calculated  from  its  definition  (4). 
Let  'll  =  ( w ,  (p ,  wn,  (p)  G  U  —  L2(Q)n  x  L2(i7)  x  Q.  Then  the  bilinear  form  in  (13)  can  be 
written  as 

b(U,  (v,  rj) )  =  -(w,  V/7  +  iujv)nh  -  (vb  iwrj  +  V  •  v  )Qh  +  {<p,  v  ■  n)dQh  +  (77,  wn)aQh. 

It  is  easy  to  check  that  in  this  case,  the  supremum  in  the  optimal  test  norm  equals 

1 1  ( v  ,V)  1 1  opt,  v  =  1 1  V?7  +  iuv  1 1  lh  +  1 1  ■ iurj  +  V  ■  v  \  \  2Qh  +  |  [v  ,  rj)  \  ^ , 


where 


=  sup 


(rj,wn)dQh  +  (^,v  -n)dQh 


(wn,fi)eQ  ||  (p) \\Q 

By  the  dehnition  of  the  norm  on  Q,  this  can  be  rewritten  as 


( z,fi)GR 


(V,z-fi)dnh  +  (ji,v  -n)dQh 

l|(A^)IU 


4.3.  Norm  Equivalence.  Now  we  turn  our  attention  to  verifying  Assumption  2.  The 
main  result  of  this  subsection  is  Theorem  4.5,  which  verifies  the  assumption.  We  be¬ 
gin  with  the  following  lemma.  Recall  that  the  generic  constant  C  is  independent  of  to 
throughout. 

Lemma  4.2:  Given  any  rj  in  L2(Q),  there  is  an  r  e  H( div,  L?)  and  a  (ft  e  Hl(Q)  such 
that 


it of  +  V  (ft  =  0, 
iuxft  +  V  •  f  —  rj, 
f  ■  n  =  ±0 


on  f2 
on  Q 
on  (9l7, 


(32a) 

(32b) 

(32c) 


(r,(ft)\\R  <  C\\rj\\n. 
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Proof.  Define  the  sesquilinear  form 

a±(</h0)  =  (V</?,  ±  uu{<p,  if)  an- 

Let  0  G  Hl(Q)  be  the  (unique)  solution  of 

a±(0,^)  =  V0  G  Hl(Q). 

Then,  set  r  by  lujf  =  — V0.  It  is  easy  to  see  that  r  G  H(div,  f2)  and  (32)  is  satisfied.  To 

prove  (33),  we  note  that  by  Assumption  3,  we  have  ||V0||^  +  w2||0||?2  <  Since 

\\iur  +V0||fife  =  0,  \\iuj(j)  +  V  -f\\nh  =  \\v\\n  and  u\\f  ||^  =  ||V0||fl,  we  immediately  have 

u2W ,  ml  =  mw^r  +  V0||2^  +  II  Ln0  +  V  •  ?\\\  +  ||f||2fi  +  ||0||2fi)  <  C^IMI^, 
which  is  (33).  □ 

The  next  lemma  is  the  analogue  of  Lemma  4.2  for  the  vector  test  variable  v  G  L2(Q)N . 

Lemma  4.3:  There  exists  an  uj\  >  0  such  that  given  any  v  in  Li2 (Q)N ,  there  is  an 
r  G  iL(div,  fi)  and  a  0  G  Hl(fi)  satisfying 

iivr  +  V 0  =  v  ,  on  i 1  (34a) 

Lu0  +  V  •  r  =  0,  on  Q  (34b) 

r  ■  n  =  ±0  on  <9i2,  (34c) 

and 

||(r,0)||R  <  C\\v\\Q  (35) 

for  all  uj  >  uj i. 

Proof.  First,  we  set  0  G  Pdl{T2)  to  be  the  (unique)  solution  of 

a±(0,^)  =  (v ,  V0)fi,  V0  G  iL1(i2).  (36) 

Then,  set  r  by  Lnr  =  — V0  +  u .  It  is  easy  to  see  that  r  G  iL(div,  1?)  and  (34)  is  satisfied. 
It  only  remains  to  prove  (35).  For  this,  we  pick  0  in  (36)  as  0  =  0  +  f  where  (  G  Hl(i 7) 
is  the  unique  solution  of  the  adjoint  problem: 

a±{<P,0  =  2cn2(</9,0)r2  V^GF1^).  (37) 

By  Assumption  3,  we  clearly  have: 

Moreover,  by  (37), 

Re (a± (0,  0  +  0)  =  Re(a±(0,  0))  +  Re(a±(0,  C)) 

=  ||V0||2r2-cn2||0||2fi  +  Re(a(0,C)) 

=  ||V0||?2  +  a;2||0||2r?. 
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Hence, 

II  V0||^  +  iv2||0||?2  =  Re(a±(0, 0  +  C))  =  (t?,V(0  +  C))fl 
<  llullfl  (l|V0||fi  +  lively) 
<C||iT||fi(||V0||2fi+a;4||^)1/2. 

\  1/2 

l  +  ^U\\lj  • 

Thus,  for  large  u,  we  obtain  ||V0||^  +  u;2||0||;^  <  Co;2||n||^.  Using  also  the  equalities 
|| iwr  +  V0|| ah  =  ||n  ||r?,  || Lu0  +  V  •  f\\nh  =  0,  and  ta||r  ||^  =  ||F  -  V0||fi,  we  obtain 

II  (r,  0)  \\r  =  II wr  +  V0||2^  +  \\h}<j>  +  V  ■  f  ||^  +  ||r  ||^  +  ||0||^ 

<  llvllfl  +  ^allu  -  V0||^+  \\</>\\n 

<  11^ II ?2  +  (ll^llfi  +  u2Ufn) 

Hence  the  estimate  of  the  lemma  follows  taking  u  large  enough.  □ 

Lemma  4.4:  There  is  an  uj\  >  0  such  that  for  any  (v,rj)  £  U  we  have 

||u||fl  +  |M|rt  <  C'||(n,77)||opty , 

for  all  ui  >  uj\. 

Proof.  For  a  given  v  and  r/,  apply  Lemmas  4.2  and  4.3  to  obtain  (r,  0)  £  R  satisfying 

iuif  +  W(f)  —  v,  on  Q  (39a) 

iu(f>  +  V  •  r  —  rj,  on  Q  (39b) 

and 

jjl  (f  ,  0)||r  <  C*  (||n  ||i7  +  |H|ri)  •  (40) 

Then 

IK7IId  +  \\v\\2n  =  far  +  V(j>,v)nh  +  (2u;0  +  V  •  r  ,  r})nh  (by  (39)), 

=  -(r  ,  iuv  +  Vrj)n  +  (V0,  v)n  +  (f ,  Vr})nh 
-  (0,  iuri  +  V  •  v)nh  +  (V  •  r  ,  rj)nh  +  (0,  V  •  v)Qh 
=  -f ,  iuv  +  Wrj)o  -  (0,  iwri  +  V  •  v  )nh 


<Cu\\v\\n  KIIV0II 


uv 


+  (0,v  -n)dQh  +  (17,  r  ■  n)doh 


(integrating  by  parts). 
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By  Cauchy-Schwarz  inequality, 


\n 


+ 


|2  ^ 


h  <  [\\VV  +  iwv\\ah  +  || iuv  +  V  •  v  U{ihj 


1/2 


ll(r,0)||* 


|  (0,  v  ■  n)anh  +  (v,  r  ■  n 
II  (r,  0)||  * 


ll(r,0)||fl. 


so  the  result  follows  from  (40).  □ 

Theorem  4.5:  Suppose  Assumption  3  holds.  Then,  there  are  positive  constants  ujq,C\,C2 
such  that  for  all  oj  >  ojq  and  for  all  (v  ,  rj)  G  H( div,  1?/,,)  x  Hl(T2h), 

Ci\\(v,v)\\v  <  II  iv  >  v)  llopt,E  <  C2\\(v,V)\\v. 

Proof. 

Upper  bound:  We  only  need  to  bound  the  jump  term  (as  the  other  terms  are  present  in 
the  Id-norm).  Since 

(/i,  v  ■  n)dQh  +  (fj,z  ■  n)dQh  =  (V//,  v)Qh  +  (n,  V  •  v  )„h  +  (z  ,  V//)^  +  (V  •  z,  rj)Qh 

=  (V/i  +  iuz,v)nh  +  (a*,^  +  V  •  u)fih 
+  (T,  V?/  +  )(2h  +  (hup,  +  V  •  z,r})(2h 

<\\(v,V)\\v\\(z,p)\\R, 


we  can  divide  by  ||(if,/i)||j?7^  0  and  take  the  supremum  on  R  to  obtain  the  upper  bound. 

Lower  bound:  We  only  need  to  estimate  the  two  terms  in  the  V -norm  that  are  not  in 
the  optimal  norm.  But  these  bounds  are  immediate  by  Lemma  4.4  provided  oj  is  large 
enough.  This  finishes  the  proof.  □ 


5.  Numerical  experiments 

In  this  section  we  present  results  of  numerical  experiments  for  four  distinct  model 
problems.  In  all  cases,  the  results  are  better  than  predicted  by  the  preceding  theory.  Let  us 
begin  by  first  describing  precisely  the  spaces  of  approximation  used  in  our  computations, 
and  the  discrete  approximation  of  the  operator  T  which  generates  the  optimal  test  space. 

We  need  to  specify  the  trial  space  (cf.  (13c)  and  (18)).  This  is  constructed  using  a  mesh 
fih  of  quadrilaterals.  Let  Q(^  denote  the  space  of  polynomials  in  one  variable  of  degree  at 
most  /,  and  let  QTm)  denote  the  space  of  polynomials  of  degree  at  most  l  and  m  in  the 
two  variables  x\  and  x-2,  resp.  We  use  it  to  define  the  sequence  of  spaces 

XV(K)  =  Qm,  YP(K )  =  Q (P*"1)  x  Q(p_1,p),  Xp_i(A0  =  Q(p-1’P-1), 

for  p  >  1.  These  sequence  of  spaces  are  to  be  interpreted  as  subspaces  of  Hl(K ), 
H(div,K),  and  L2(K),  resp.,  where  K  =  (0,  l)2.  The  spaces  Xp(K),  Yp(K ),  Xp_i(K) 
on  a  general  K  G  L2h  are  the  analogous  spaces  of  shape  functions  on  the  physical  element 
defined  through  the  usual  pullback  (depending  on  the  subspace  interpretation)  to  the  refer¬ 
ence  element.  Similarly,  let  Q ^(E)  on  a  (possibly  curved)  mesh  edge  E  denote  the  image  of 
Q W  from  (0, 1)  under  the  standard  map.  The  interior  variables  are  approximated  in  Sh  = 
{r  :  t\k  ^  XP(K)N}  and  Wh  =  {v  :  v\k  £  XP(K)},  resp.,  while  the  numerical  fluxes  and 
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traces  lie  in  Qh  =  {(zUth,  hh)  G  R  '■  4,a|e  and  jlh\E  are  in  Q(p+y)(E)  for  all  mesh  edges  E, 
and  zU:h  —  fih  =  0  on  d i?},  i.e., 

Uh  =  ShxWhx  Qh. 

Note  that  the  numerical  traces  (in  the  second  component  of  Qh)  are  continuous  at  edges 
that  meet  at  a  vertex. 

As  noted  in  §  2.3,  the  practical  application  of  DPG  requires  us  to  approximate  the 
operator  T  by  a  discrete  version  T  which  maps  into  a  finite  dimensional  enriched  test 
space  V  C  V.  In  our  experiments,  V  is  constructed  by  considering  the  local  polynomial 
order  p  of  the  element  K  and  a  global  parameter  A p.  We  compute  the  practical  test  space 
using  (10),  with  V  set  element  by  element  by 

V\K  =  Yp+Ap(K)  x  Xp+Ap(K)  C  H(div,  K)  x  H\K). 

In  all  our  experiments,  we  fix  A p  =  2.  This  choice  is  dictated  by  our  previous  computa¬ 
tional  experience  [6,  8,  9],  whereby  it  was  clear  that  using  a  higher  A p  did  not  result  in 
any  error  improvements,  while  using  a  lower  A p  could  result  in  a  non-injective  T. 

5.1.  Model  problem  A.  This  is  problem  (11),  for  which  we  provided  a  theoretical  anal¬ 
ysis.  Our  domain  Q  is  the  square  (0,  l)2.  The  right  hand  sides  /  and  g  of  (lib)  and  (11c), 
respectively,  are  set  so  that  the  exact  solution  is  a  plane  wave  solution  (propagating  in 
the  ^-direction) ,  namely 

0(£)  =  e-Mzicos^sin^  =  0(f)  (cos  6,  sind), 

and  correspondingly, 

/  =  0,  g(x)  =  (f>(x)[rii(x)  cos 9  +  n2(x)  sin#  —  1] 

where  the  outward  unit  normal  vector  is  n  =  (711,712). 

We  employ  a  grid  of  square,  bilinear  elements  (with  h  denoting  the  length  of  their  side) 
to  discretize  the  problem,  i.e.,  p  —  1  in  the  above  definition  of  Uh-  Our  mesh  is  quite 
coarse,  with  uh  =  n/2,  equivalent  to  four  elements  per  wavelength  when  the  direction 
is  aligned  with  the  mesh.  This  is  still  enough  to  reasonably  represent  the  wave,  as  the 
L2  best  approximation  error  (BAE)  will  vary  from  about  6.5%  to  9%,  depending  on  the 
direction  6  (the  case  of  diagonal  propagation  6  =  7t/4,  being  the  best  case).  We  use  BAE 
to  denote  the  best  approximation  error  in  both  u  and  0,  i.e., 

(BAE)2  =  inf  || u  -  wh \\2n  +  inf  ||0  -  <ph \\2ni 
whesh  ‘Ph&Wh 

while  DE  denotes  the  discretization  error  defined  by 

(DE)2  =  || u  -  uh\\2n  +  ||0  -  0fc||^. 

We  will  report  the  ratio  DE/BAE.  Due  to  Theorem  3.1,  we  expect  this  ratio  to  be  bounded 
by  a  mesh- independent  constant.  We  are  interested  in  how  this  ratio  changes  with  u.  Since 
BAE  represents  the  best  any  method  can  do  using  the  given  spaces,  a  study  of  how  the 
ratio  DE/BAE  changes  with  respect  to  u  can  give  us  an  idea  of  how  pollution  effects 
influence  discretization  errors. 

We  compare  the  performance  of  the  DPG  method  to  two  other  methods  readily  available 
in  our  software  package,  namely  the  standard  FEM,  and  the  FEM  using  the  new  and 
interesting  quadrature  rules  of  [1].  Note  that  when  considering  these  other  methods,  BAE 
is  to  be  interpreted  as  the  L 2  best  approximation  error  of  the  space  used  by  the  other 
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DE/BAE  with  four  bilinear  elements  per  wave  (0=0)  DE/BAE  with  four  bilinear  elements  per  wave  (0=jt/4) 


FIGURE  1.  Model  A:  The  ratio  of  discretization  errors  to  best  approxima¬ 
tion  erros  in  L2(l?)-norrn  for  two  plane  wave  directions. 

method,  i.e. ,  BAE=  inf^g^i^nw^  ||^> — V?*.|| r?,  and  similarly,  DE  denotes  the  discretization 
error  of  the  other  method.  Note  also  that  we  use  a  very  modest  number  of  mesh  points 
per  wavelength  (e.g.,  in  comparison  to  results  from  other  least  square  methods  in  the 
literature  [20]). 

The  results  are  in  Figure  1  for  two  values  of  9.  We  observe  that  in  both  cases  the 
quality  of  standard  finite  element  approximations  quickly  deteriorates  as  we  increase  the 
wavenumber.  The  deterioration  exists,  but  is  delayed,  when  the  method  of  [1]  is  used  - 
see  the  curve  labeled  “AW  Quadrature” .  The  DPG  solutions  however  are  still  very  close 
to  the  best  approximations. 

5.2.  Model  problem  B.  This  problem  is  similar  to  Model  Problem  A.  We  keep  all 
parameters  the  same  as  before,  but  use  slightly  different  boundary  conditions.  These 
boundary  conditions  are  typically  used  to  demonstrate  the  accumulation  of  phase  errors 
in  the  direction  of  wave  propagation  for  standard  Galerkin  methods.  The  boundary  value 
problem  is: 


ILUU  +  V0  =  0, 

on  i? 

(41a) 

lL0(j)  +  V  ■  U  =  /, 

on  i? 

(41b) 

u  ■  n  =  un, 

on  Ti  U  r4 

(41c) 

u  ■  n  -  0  =  g, 

on  r2  ur3, 

(41d) 

i.e.,  we  prescribe  the  normal  “velocity” 

on  the  lower  and  left  edges  of  the  domain 

(riu 

r4)  and  maintain  the  previous  Robin  boundary  conditions  at  the  upper  and  right  edges 
(r2  U  r3).  We  choose  /,  un,  and  g  so  that  the  exact  solution  is  the  same  as  in  Model 
Problem  A. 

The  results  are  depicted  in  Figure  2.  In  Figure  2(a),  we  observe  a  behavior  similar  to 
Figure  1.  In  Figure  2(b),  we  see  that  the  relative  L2  error  percentage  remains  more  or  less 
constant  (about  8%)  for  the  DPG  method,  while  it  increases  with  the  number  of  degrees 
of  freedom  to  about  140%  for  the  standard  method.  (Note  that,  as  before,  the  number  of 
degrees  of  freedom  is  tied  to  the  wavenumber  through  uh  =  n/2.) 
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DE/BAE  with  four  bilinear  elements  per  wave  (0=ji/4) 


Wavenumber  co  (on  log  scale) 


(b)  Relative  L 2  error  (in  %)  for  the  three  methods 


1 1.00 
0.750 
0.500 
0.250 


(c)  Left:  Plot  of  the  error  for  the  CG  method  showing  accumulation  of  phase  errors  (at 
the  northeast  corner)  in  the  direction  of  wave  propagation.  Right:  Plot  of  the  error  for  the 
DPG  method  showing  no  accumulation  of  phase  errors  in  the  propagation  direction. 

Figure  2.  Results  for  Model  B  (all  results  are  for  6  =  7t/4) 


5.3.  Model  problem  C.  In  this  problem,  the  exact  solution  consists  of  the  cylindrical 
wave 

0(f)  =  H^\u\x\), 

(2) 

where  H q  ;  is  the  zero-order  Hankel  function  of  the  second  kind.  The  domain  consists 
of  the  square  ( — 1,  l)2  with  a  circular  exclusion  of  radius  a  =  0.1  in  the  center,  i.e. 
ft  =  (—1,  l)2  \  {x  :  \x |  <  a}.  We  denote  the  boundary  of  the  circle  by  Ta.  The  equations 
of  the  boundary  value  problem  are 


icon  +  V0  =  0, 

on  ft 

(42a) 

iuxj)  +  V  •  u  =0, 

on  ft 

(42b) 

u  ■  n  —  un, 

on  Ta 

(42c) 

u  •  n  -  0  =  g, 

on  dft  \  Ta. 

(42d) 
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(a)  Mesh  when  w  =  I67 r.  Note  the  use  of 
curved  elements. 


10  25  50  100  200  400 

Wavenumber  co  (on  log  scale) 


(c)  The  DPG  method,  at  its  lowest  order  (p  =  1), 
compares  favorably  to  other  higher  order  ( p  =  2) 
methods. 


DE/BAE  when  to  h  =  jt/2 


10'  102  103 


Wavenumber  to  (on  log  scale) 

(b)  Ratio  of  the  discretization  error  to  the  best  ap¬ 
proximation  error  (in  L2-norm)  for  three  methods. 


DE/BAE  for  DPG  method  when  to  h  =  jt/2  and  p=1 


(d)  A  comparison  of  the  performance  of  DPG 
method  for  this  model  and  the  previous  two  model 
problems  (with  zoomed  ordinate  scale). 


Figure  3.  Results  for  Model  C. 


We  discretize  using  a  grid  of  square  elements,  with  the  exception  of  a  block  of  elements 
which  are  deformed  by  geometry  mappings  generated  through  transfinite  interpolation  - 
see  Figure  3(a). 

The  results  are  similar  to  the  previous  two  model  problems.  The  DPG  error  closely 
follows  the  L 2  best  approximation  error,  as  evidenced  by  the  ratio  plotted  in  Figure  3(b), 
which  remains  close  to  the  optimal  value  of  1.  These  results  are  for  p  =  1. 

We  also  compared  these  results  with  the  p  =  2  case  of  the  standard  FEM.  It  is  well 
known  that  increasing  the  order  improves  pollution  errors  for  the  standard  method.  This 
is  indeed  the  case,  as  seen  in  Figure  3(c),  for  both  the  standard  FEM  and  its  modification 
of  [1],  However,  even  with  this  improved  performance,  neither  of  these  methods  compared 
favorably  to  the  lowest  order  DPG  method.  While  the  DPG  error  remained  under  9% 
(in  the  L 2  norm)  for  the  range  of  wave  numbers  considered,  the  L2  errors  of  the  other 
methods  eventually  increased  beyond  a  100%. 


DPG  METHOD  FOR  WAVES 


21 


Finally,  in  Figure  3(d),  to  get  an  idea  of  the  relative  difficulty  of  the  model  problems 
we  considered  so  far,  we  compared  the  performance  of  the  DPG  method  for  Models  A, 
B,  and  C.  The  DE/BAE  ratio  remains  close  to  the  optimal  value  of  one  for  all  the  three 
cases.  Model  A,  run  with  the  propagation  direction  9  —  0  in  this  plot,  seems  to  show  the 
largest  increase  in  DE/BAE  as  the  wavenumber  increases.  (We  have  run  this  model  to 
the  limit  of  our  computational  resources.)  Although  this  increase  is  a  small  fraction  of 
the  increase  the  other  methods  suffer  from,  the  fact  that  there  is  a  slight  increase  seems 
to  indicate  that  the  DPG  method  may  also  suffer  from  pollution  errors  for  large  enough 
wave  numbers.  The  current  data  however  is  insufficient  to  make  a  definitive  conclusion. 

5.4.  Model  problem  D:  Pekeris  waveguide.  Finally,  we  consider  a  more  realistic 
example  of  wave  propagation.  The  Pekeris  waveguide  (see  Figure  4(a))  is  a  canonical 
example  of  a  shallow  water  waveguide.  This  model  consists  of  a  water  layer  above  a 
sediment  layer.  A  point  source  within  the  water  column  at  depth  zs  generates  time- 
harmonic  pulses  which  propagate  into  the  water  and  sediment  layers.  The  sediment  layer 
is  represented  as  an  acoustic  medium  with  higher  density  and  sound  speed.  The  change  in 
acoustic  properties  occurs  at  a  depth  H.  At  the  surface  Ti  of  the  waveguide,  a  pressure- 
release  boundary  condition  is  prescribed  (i.e.  0  =  0).  The  speed  of  sound  in  water  is 
taken  to  be  c  =  1500  meters  per  second.  We  set  L  =  1500  m  to  be  our  length  scale  in 
non-dimensionalization  of  the  problem.  An  additional  scaling  is  applied  to  the  ambient 
density  so  that  p0  =  1  within  the  water  column.  The  result  of  this  scaling  is  such  that 
u  and  0  are  of  the  same  order  of  magnitude.  The  full  set  of  problem  parameters  after 
non-dimensionalization  is  as  follows: 


Cl  =  1 

speed  of  sound  in  water 

c2  =  1.2 

speed  of  sound  in  sediment 

pi  =  1 

ambient  density  of  water 

p-2  =  1-8 

ambient  density  of  sediment 

H  =  1/15 

depth  of  water  column 

=  36/1500 

depth  of  point  source. 

The  original  Pekeris  problem  is  posed 

on  the  unbounded  half  space  {(x,z) 

C  M2  : 

z  <  0}  with  Sommerfeld  radiation  conditions.  This  can  be  numerically  solved  using 
PML  or  other  techniques  for  truncating  infinite  domain  problems.  However,  since  these 

truncation  techniques  are  not  the  subject  of  the  present  study,  we  construct  a 

simpler 

model  problem  by  truncating  to  the  domain  D  =  [0,  200/L]  x  [0,  H/L\  and  simply  imposing 
non-homogeneous  Robin  boundary  conditions  with  contrived  data  obtained  from  the  exact 

solution: 

ioo  _ 

0  +  V  -u  =  dZs , 

in  12 

(43a) 

PoC2 

icopou  +  V0  =  0, 

in  12 

(43b) 

o' 

-e- 

on  r, 

(43c) 

U  •  n  -  0  =  g, 

on  r2 

(43d) 

u  ■  n  —  un, 

on  r3 

(43e) 
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rx 


X 


Water 

O  Point  source 

Cl  =  l,Pl  =  1 

Sediment 

C2  =  1.2,  P2  =  1.8 

r3 

(a)  Diagram  of  the  Pekeris  waveguide. 


(b)  The  exact  solution  for  lj  =  2n  x  400. 


(c)  Error  from  standard  biquadratic  FEM  using  (d)  Solution  from  standard  biquadratic  FEM  using 
about  4  elements  per  wavelength  (ut  =  2tt  x  400).  about  4  elements  per  wavelength  (to  =  27t  x  400). 


(e)  Error  from  DPG  method  using  about  4  bilinear  (f)  Solution  from  DPG  method  with  4  bilinear  ( p  = 
elements  per  wavelength  (lo  =  2ir  x  400).  1)  elements  per  wavelength  (uj  =  2n  x  400). 


Figure  4.  Model  D:  Pekeris  waveguide.  Only  the  real  parts  of  the  scalar 
component  are  shown.  The  color  scale  for  Figures  4(b),  4(d),  and  4(f)  is 
[—10  •  ■  •  10],  while  for  Figures  4(c)  and  4(e),  it  is  [0  ■  •  •  5]. 


An  analytic  expression  for  the  pressure  0  within  the  water  column,  which  can  be  derived 
through  Fourier  transformation  and  complex-contour  integration  (see,  e.g.  [19]),  is  given 
by 


(p(x,z) 


_2?  r  7i72  sin(7iza)  sin^z)  , 

AAFM  7 fg  cos2(7ltf)  +  7lg  sin2(7ltf)  '  71 

N 

+  ^  um(zs)um(z)e~l0rnX , 

PKzs)  ^ 


(44) 
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where 

'  kx  =  u  /  Ci, 

k2  =  u/c2, 

P  =  Vki  7i  >  0, 

J2  =  \Jk\  —  P2  >  0, 

and  the  N  modal  values  =  \Jk\  —  y2m  =  VM  —  7 2m  >  0  are  determined  by  the  N 
solutions  of  the  characteristic  equation 

,  /  tt\ 

tan(7i  mH)  —  - . 

72mPl 

A  plot  of  the  exact  solution  is  shown  in  Figure  4(b). 

Plots  of  the  errors  and  computed  solutions  are  in  Figure  4.  Since  the  results  from 
the  lowest  order  (bilinear)  case  of  the  FEM  are  poor  (not  displayed),  we  only  compare 
the  bilinear  DPG  solution  to  the  biquadratic  FEM.  A  close  observation  of  the  solution 
computed  using  the  standard  FEM  in  Figure  4(d)  shows  a  small  phase  lag  when  compared 
to  the  exact  solution  in  Figure  4(b).  Since  this  may  be  hard  to  judge  from  Figure  4(b), 
we  also  include  visualization  of  the  errors  in  Figures  4(c)  and  4(e).  Clearly,  away  from  the 
source,  the  error  is  large  for  the  standard  method,  while  for  the  DPG  method,  the  error 
remains  more  or  less  the  same  throughout  the  domain.  Note  also  that  small  phase  errors 
can  lead  to  large  L 2  errors.  The  DPG  solution,  even  in  the  lowest  order  case,  shown  in 
Figure  4(f),  is  a  far  better  approximation  (visually  almost  identical  to  the  exact  solution). 

6.  Conclusion 

6.1.  Summary.  We  presented  a  new  DPG  method  for  acoustic  time  harmonic  wave  prop¬ 
agation.  Although  this  method  has  more  unknowns  than  other  standard  methods,  and 
although  it  does  not  have  conservation  properties  in  its  current  form,  we  think  it  is  an 
interesting  alternative  because  it  exhibits  remarkably  small  phase  errors  in  all  attempted 
numerical  experiments.  While  many  standard  methods  show  comparable  performance  for 
low  to  moderate  wave  numbers,  for  large  wave  numbers,  the  new  method  is  highly  com¬ 
petitive.  Our  analysis  using  the  known  regularity  and  stability  results  for  the  Helmholtz 
equation  leads  to  a  proof  of  error  estimates,  which  however  does  not  explain  the  low  phase 
errors. 

6.2.  The  analysis  in  hindsight.  One  may  observe,  as  we  did  in  hindsight,  that  the 
analysis  we  performed  in  Section  4  has  elements  that  can  be  generalized  to  apply  for 
various  problems  beyond  wave  propagation.  To  briefly  remark  on  a  way  to  generalize, 
consider  any  abstract  problem  Au  =  /,  with  a  linear  operator  A  :  D(A)  i— >•  L2(i7),  where 
the  domain  of  A,  D(A),  incorporates  any  boundary  conditions  on  u,  and  is  equipped  with 
the  graph  norm  ||w|||,^x  =  ||u||^  +  ||Au||^.  Formally  introducing  an  L2  adjoint  operator 
A*  by  (Au,  v )  12  =  (u,  A*v) nh  +  ((u,  v))gnh,  where  we  have  lumped  all  element  interface  and 
boundary  terms  into  ((•,  -))aah,  we  can  P°se  an  abstract  ultraweak  formulation  as  follows: 
Find  u  in  L2(i?)  satisfying  (u,A*v)ah  +  (( u,v))aah  =  (f,v)nh  for  all  v  in  the  space  V  with 
norm  defined  by  H2V  =  IMI^  +  P*^-  If  we  set  A  to  the  Helmholtz  (Erst  order)  wave 
operator,  we  End  that  this  is  exactly  the  formulation  on  which  the  DPG  method  of  this 
paper  is  based,  cf.  (13).  Note  that  in  the  above  generalization,  u  is  sought  in  a  space  Q, 
which  is  normed  by  a  straightforward  generalization  of  the  quotient  norms  in  §  2.4.  We 
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can  then  view  the  entire  analysis  of  Section  4  as  aimed  at  proving  the  equivalence  of  the 
||  •  ||y-norm  with  the  “optimal  norm” 


v 


2 

opt,V 


\\A*v\\lh  +  |  MIL,’  where 


sup 

u£D(A) 


{{u,v))anh 

IMb(A) 


cf.  §4.2.  The  inequality  ||w||opty  <  c\\v\\  v  can  be  proved,  even  in  the  general  context, 
exactly  as  in  proof  of  the  upper  bound  in  Theorem  4.5.  The  gist  of  the  argument  to 
prove  the  reverse  inequality  can  be  abstracted  from  Section  4,  under  the  assumption  that 
IM|j7  <  Co 1 1 Aw H/2.  (This  assumption  follows  from  Lemmas  4.2  and  4.3  in  our  particular 
case.)  Then,  given  any  v  E  V,  considering  a  u  that  solves  Au  =  v,  we  have 

INI??  =  (Au,v)n  =  (u,A*v)nh  +  ((u,v))ds 4 

<  ||w||d|| A*n||/2h  +  ||«||u(A)|[u]|aflh 

^  (^olMli?  +  (Co  +  i)IMIi?)  (M*vII?4  +  |M la/2h)  ’ 

which  proves  that  ||w||y  <  C  ||w||  t  v,  thus  completing  the  proof  of  the  norm  equivalence. 
In  hindsight,  we  understand  that  this  is  the  essence  of  the  analysis  not  only  in  this  paper, 
but  also  in  [6]. 


6.3.  Future  directions.  Our  future  studies  are  focused  on  sharpening  our  theoretical 
tools  and  formalizing  and  improving  the  informally  stated  abstract  generalization  above 
(§  6.2),  all  with  the  aim  of  proving  or  disproving  that  the  errors  are  pollution-free.  Another 
direction  of  research  we  are  exploring  exploits  the  Hcrmitian  positive  definiteness  of  the 
linear  system,  together  with  the  good  phase  errors,  to  design  efficient  iterative  solvers. 
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